URD

Load libraries.

library(Seurat)

Registered S3 method overwritten by 'data.table':
  method           from
  print.data.table     
Registered S3 method overwritten by 'htmlwidgets':
  method           from         
  print.htmlwidget tools:rstudio
Attaching SeuratObject

library(URD)

Loading required package: ggplot2
Loading required package: Matrix
Registered S3 method overwritten by 'gplots':
  method         from 
  reorder.factor gdata

library(tidyverse)

Registered S3 methods overwritten by 'dbplyr':
  method         from
  print.tbl_lazy     
  print.tbl_sql      
── Attaching packages ────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse 1.3.2 ──✔ tibble  3.1.8      ✔ dplyr   1.0.10
✔ tidyr   1.2.1      ✔ stringr 1.5.0 
✔ readr   2.1.3      ✔ forcats 0.5.2 
✔ purrr   1.0.1      ── Conflicts ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
✖ tidyr::expand() masks Matrix::expand()
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
✖ tidyr::pack()   masks Matrix::pack()
✖ tidyr::unpack() masks Matrix::unpack()

Load merged eb, lb, ehf.

load("/Volumes/Bud_SSD/20230122-urd/merged-seurat.Robj")

Grab count data and meta data.

count.merged <- as.matrix(merged.seurat@assays$RNA@counts)

Warning: sparse->dense coercion: allocating vector of size 6.2 GiB

meta.merged <- as.data.frame(merged.seurat@meta.data)

rm(merged.seurat)

Create urd.

merged.urd <- createURD(count.data = count.merged, meta = meta.merged, min.cells = 3, min.counts = 0)

2023-01-23 07:20:42: Filtering cells by number of genes.
2023-01-23 07:20:59: Filtering genes by number of cells.
2023-01-23 07:21:23: Filtering genes by number of counts across entire data.
2023-01-23 07:21:46: Filtering genes by maximum observed expression.
2023-01-23 07:22:11: Creating URD object.
2023-01-23 07:22:15: Determining normalization factors.
Warning: sparse->dense coercion: allocating vector of size 6.0 GiB2023-01-23 07:22:28: Normalizing and log-transforming the data.
as(<dgeMatrix>, "dgCMatrix") is deprecated since Matrix 1.5-0; do as(., "CsparseMatrix") instead
2023-01-23 07:23:12: Finishing setup of the URD object.
Warning: sparse->dense coercion: allocating vector of size 6.0 GiB2023-01-23 07:23:26: All done.

Add stage.

merged.urd@group.ids$stage <- as.character(merged.urd@meta[rownames(merged.urd@group.ids),"stage"])

Add embryo.

merged.urd@group.ids$embryo <- as.character(merged.urd@meta[rownames(merged.urd@group.ids),"embryo"])

#embryo <- sort(unique(merged.urd@group.ids$embryo))
embryo <- c("EB1", "EB2", "AVJ4-56", "AYQ8-11", "AYR4-57", "BBH6-17", "AXL7-21", "EHF1", "EHF2")

Calculate variable genes. Use only Flox embryos.

var.by.embryo <- lapply(c(1, 2, 3, 4, 5, 6, 7, 8, 9), function(n) {
  findVariableGenes(merged.urd, cells.fit=cellsInCluster(merged.urd, "embryo", embryo[n]), set.object.var.genes=F, diffCV.cutoff=0.3, mean.min=.005, mean.max=100, main.use=paste0("Embryo ", embryo[n]), do.plot=T)
})

Add variable genes.

var.genes <- sort(unique(unlist(var.by.embryo)))
merged.urd@var.genes <- var.genes

Calculate PCA.

merged.urd <- calcPCA(merged.urd, mp.factor = 2)

[1] "2023-01-23 07:32:26: Centering and scaling data."
[1] "2023-01-23 07:32:34: Removing genes with no variation."
[1] "2023-01-23 07:32:36: Calculating PCA."
[1] "2023-01-23 07:37:28: Estimating significant PCs."
[1] "Marchenko-Pastur eigenvalue null upper bound: 1.40724420496983"
[1] "31 PCs have eigenvalues larger than 2 times null upper bound."
[1] "Storing 62 PCs."

pcSDPlot(merged.urd)

set.seed(19)
merged.urd <- calcTsne(object = merged.urd)
plotDim(merged.urd, "stage", plot.title = "tSNE: Stage")

Warning: The `<scale>` argument of `guides()` cannot be `FALSE`. Use "none" instead as of ggplot2 3.3.4.

plotDim(merged.urd, "cluster.whole", plot.title = "tSNE: Cluster")

Save.

save(merged.urd, file = "merged-tsne-urd.Robj")

Load.

load("merged-tsne-urd.Robj")

Calculate diffusion map.

merged.urd <- calcDM(merged.urd, knn = sqrt(dim(merged.urd@meta)[1]), sigma=NULL)

[1] "destiny determined an optimal global sigma of 41.63"

Warning: You have 1603 genes. Consider passing e.g. n_pcs = 50 to speed up computation.as(<dsCMatrix>, "dgTMatrix") is deprecated since Matrix 1.5-0; do as(as(., "generalMatrix"), "TsparseMatrix") instead

plotDimArray(merged.urd, reduction.use = "dm", dims.to.plot = 1:8, outer.title = "Diffusion Map (Sigma NULL, 214 NNs): Stage", label="stage", plot.title="", legend=F)

Warning: The `<scale>` argument of `guides()` cannot be `FALSE`. Use "none" instead as of ggplot2 3.3.4.

plotDim(merged.urd, "stage", transitions.plot = 10000, plot.title="Developmental stage (with transitions)")

plotDim(merged.urd, "cluster.whole", transitions.plot = 10000, plot.title="Developmental stage (with transitions)")

Add cluster and stage.

merged.urd@meta <- merged.urd@meta %>% unite("cluster.stage", cluster.whole:stage, remove = FALSE)
merged.urd@group.ids$cluster.stage <- as.character(merged.urd@meta[rownames(merged.urd@group.ids),"cluster.stage"])

Calculate pseudotime.

pseudotimePlotStabilityOverall(merged.urd)

plotDim(merged.urd, "pseudotime")


plotDists(merged.urd, "pseudotime", "stage", plot.title="Pseudotime by stage")

plotDists(merged.urd, "pseudotime", "embryo", plot.title="Pseudotime by embryo")

plotDists(merged.urd, "pseudotime", "genotype", plot.title="Pseudotime by genotype")

Save.

save(merged.urd, file = "merged-pseudotime-urd.Robj")

Extract pseudotime.

pseudo.data <- data.frame(rownames(merged.urd@pseudotime), merged.urd@meta$stage, merged.urd@meta$embryo, merged.urd@meta$genotype, merged.urd@meta$cluster.whole, merged.urd@meta$germ.whole, merged.urd@pseudotime$pseudotime)

save(pseudo.data, file = "pseudotime-urd.Robj")